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The genetic diversity of a species is shaped by its recent evolutionary history and can be used 
to infer demographic events or selective sweeps. Most inference methods are based on the null 
hypothesis that natural selection is a weak evolutionary force. However, many species, particularly 
pathogens, are under continuous pressure to adapt in response to changing environments. A 
statistical framework for inference from diversity data of such populations is currently lacking. 
Toward this goal, we explore the properties of genealogies that emerge from models of continual 
adaptation. We show that lineages trace back to a small pool of highly fit ancestors, in which 
simultaneous coalescence of more than two lineages frequently occurs. While such multiple mergers 
are unlikely under the neutral coalescent, they create a unique genetic footprint in adapting 
populations. The site frequency spectrum of derived neutral alleles, for example, is non-monotonic 
and has a peak at high frequencies, whereas Tajima's D becomes more and more negative with 
increasing sample size. Since multiple merger coalescents emerge in various evolutionary scenarios 
characterized by sustained selection pressures, we argue that they should be considered as a null- 
model for adapting populations. 



Evolutionary change often takes longer than the win- 
dow of observation. A sequence sample therefore rep- 
resents a static snapshot from which we want to learn 
about a dynamic evolutionary process. The predominant 
framework to analyze population genetic data and infer 
demographic history is the neutral coalescent. Within 
the neutral coalescent, all individuals are equivalent, i.e., 
there are no fitness differences. The statistical proper- 
ties of the genealogies that arise in this simple popu- 
lation genetic model can be computed exactly (King- 



pathogens like influenza that continuously evade human 



man 



1982 Wakeley 2008), facilitating comparison to 



data. One central prediction of the neutral coalescent 
is that the genetic diversity of a population is propor- 
tional to its size. This prediction, however, is at odds 
with the observed weak correlation between genetic di- 
versity and population size, a paradox often remedied by 
the definition of an effective population size via the ge- 
netic diversity (Charlesworth 20091. But the model is 



also too simple in other aspects and therefore has been 
generalized to account for historic changes in population 
size, mutation rates, geographical structure and the ef- 



fects of purifying selection 


( Barton and Etheridge 2004 


Charlesworth et al. 1993 Nordborg 1997 O'fallon et al. 


2010 


Walczak et al. 


2012) 


. Positive selection, however 



has proved difficult to incorporate, and progress has been 



and Schweinsberg 


Neuhauser 


1997). 



In many populations, particularly large microbial pop- 
ulations, selection is neither rare nor weak. Instead, 
these populations are under sustained pressure to adapt 
to changing environments. Prominent examples include 



immune responses ( Smith et al. 2004 ) or HIV, which es- 



tablishes a chronic infection despite heavy immune pre- 
dation ( [Lemey et al. 2006). The genealogical trees recon- 
structed from sequence samples often suggest substantial 
departure from neutrality: see ( [Bedford et al 



2011) 



2010) 



for 
for 



examples from viral evolution or (Seger et al. 
cukaryotic examples. The influenza tree shown in Fig. [I] 
for instance, is incompatible with a neutral genealogy, 
since there are parts where many lineages merge in a very 
brief period, and the tree often branches extremely un- 
evenly with very few individuals on the right and many on 
the left. These two observations represent fundamental 
deviations from the standard neutral model, even when 
a varying population size is allowed. 

To analyze and interpret genealogies of populations un- 
der sustained directional selection, an alternative simple 
null model would be extremely useful. The features of ge- 
nealogies discussed above are in fact common to a class 
of non-Kingman like coalescence models, which have re- 
ceived considerable attention in the mathematical coales- 



cent literature (Berestycki 2009 Pitman 19991 



A special case of these non-Kingman coalescent 
models, the Bolthausen-Sznitman coalescent (BSC) 
( Bolthausen and Sznitman[|1998[ ), has been shown to de- 
scribe the genealogies in models where a population ex- 
pands into uninhabited territory ( jBrunet et «/.||2007[ ). In 
this case, individuals at the front of the population enjoy 
a growth advantage. As a result, the population develops 
moving front ( |Fisher[ |1937| |Kolmo gorov et al.l 19371) 



In a remarkable tour de force, Brunet et al. have solved a 



special case of such models exactly (Brunet et al. 2007) 



2 



skewed branching 
multiple mergers 




ber of offspring is one. Similar models have been used 



FIG. 1 Panel A shows a maximum-likelihood tree of influena 
sequences (HA segment) sampled in Asia in 2009 (subtype 
H3N2) produced using Fasttree ( |Price et ai]\2009\ . Panel B 



shows a tree drawn from a simulation of the model. 



and provided a phenomcnological theory making the case 
that genealogies obey the same statistics for all mo dels 
of this kind (FKPP); see ( |Brunet and Derrida[ [2012] ) for 
a recent review. 



We will argue in this article that the Bolthausen- 
Sznitman coalescent (BSC) emerges generically in models 
of rapidly adapting populations in a similar way as it de- 
scribes genealogies in traveling waves of FKPP type. We 
perform extensive computer simulations and investigate 
the distribution of heterozygosity in the population, the 
average time to the most recent common ancestor, and 
the site frequency spectrum (SFS). Most notably, the site 
frequency spectrum is non-monotonic with a large num- 
ber of high frequency derived alleles, as is often found in 
samples of adapting populations. We then present heuris- 
tic arguments that the underlying genealogical process is 
approximately the BSC. In the discussion, we outline the 
basic features of the BSC and discuss its applicability to 
wider classes of models. The BSC may thus be consid- 
ered as a null-model for populations in which turnover is 
driven by selection, i.e., where genetic drift is negligible 
compared to genetic draft or selective interference. 



I. THE MODEL 



by a number of authors (Desai and Fisher 


2007 


Hal- 


latschek[ 2011 


Neher et a . 2010 Park and Krug| |2007| 


Rouzine et al. 


2003 


Tsimring et al.\ 


19961. In our model, 



the fitness of offspring is given by their parents' value in- 
cremented by a Gaussian random number with variance 
s 2 . The latter process mimics mutations that change 
the fitness by small and frequent increments. Quite gen- 
erally, the competition between dispersal of individuals 
along the fitness axis (mutation) and selection results in 
a population that behaves as a traveling pulse along the 
fitness axis; see Fig. [2] Absolute fitness itself is of course 
not increasing indefinitely, but increasing fitness is offset 
by environmental deterioration. 

Most of the above-mentioned investigations focused on 
the rate of adaptation, which is rarely observable. In- 
stead, we focus on genealogies and their relation to ob- 
served genetic diversity. Our model can be easily im- 
plemented as a computer program that keeps track of 
the complete genealogy of the population. It also has 
a straightforward continuous time limit described by a 
stochastic differential equation for the concentration pro- 
file c(x, i) 

dtc(x, t) = Dd 2 c(x, t) + (x — X)c(x, t) + y/ c(x, t)i](x. t) , 

where the last term accounts for the stochastic nature of 
reproduction. The diffusion constant D = (is 2 / 2 (the 
mutation rate fj, = 1 in the discrete model discussed 
above). The term A is given by the product of the popu- 
lation averaged growth rate and N/Nq. The latter keeps 
the population size approximately constant at Nq. 

In this model, large populations attain a steady fitness 
distribution of roughly Gaussian shape, which translates 
with a velocity v rj Ds (24 log N)~s towa r ds higher fitness, 



where N = NDa (Cohen et al. 



2005 



Tsimring et al. 



1996]). The distribution and its landmarks are sketchec. 
in Fig. [2] It is convenient to measure fitness relative to 
the population mean, x — x — x; we will drop the tilde in 
the following. The fittest individuals of the population 
reside roughly x c « v 2 /AD above the population mean. 
Computer programs and analysis scripts are available on 
request. 



The evolutionary dynamics of a population are mainly 
determined by the distribution of fitness in the popula- 
tion, where fitness is the expected number of offspring 
an individual leaves behind. In general, fitness depends 
on many traits, which are affected by mutations. In a 
rapidly changing environment, populations are far from 
any fitness optimum, with many mutations available that 
increase fitness (and even more that decrease fitness). 

To model such scenarios, we consider a collection of 
N individuals that are characterized by a (log-)fitness 
x, which determines their average reproductive success. 
Specifically, the number offspring of an individual is Pois- 
son distributed with mean exp(x— x), where x is the mean 
fitness, which assures that the population averaged num- 



II. RESULTS 

We first present simulation results of our model and 
contrast the patterns of genetic diversity of continuously 
adapting populations with neutral expectations. Below, 
we will analyze our model mathematically and show that 
the striking differences result from the exponential am- 
plification of individual lineages by selection. 

Fig. [Tj3 shows the genealogy of a sample drawn from 
our simulations. This tree is characterized by long termi- 
nal branches, very uneven branching of the deep nodes, 
and bursts of rapid merging of lineages. We will investi- 
gate these features in the following. 
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FIG. 2 Ancestral lineages in evolving populations. The fig- 
ure shows the fitness distribution of the population, translat- 
ing towards higher fitness with velocity v, at two time points. 
Randomly sampled individuals (green, blue, and violet dots in 
the later population) tend to come from the center of the dis- 
tribution, while ancestors tend to be among the fittest in the 
population. The ancestral lineages wiggle due to mutations 
that randomly perturb their fitness. Simultaneously, lineages 
move towards the high fitness edge, where they are likely to 
meet and coalesce. The fittest individuals are typically at 
x c ~ v 2 /AD above the mean fitness. 



A. Distribution of heterozygosity and pair coalescence 




FIG. 3 The distribution of pair coalescence times (propor- 
tional to heterozygosity) in a model of rapidly adapting pop- 
ulations. Note that the time axis is rescaled by The 
fact that all curves for different N and D collapse onto a sin- 
gle master curve demonstrates that oc Ds logs N is the 
time scale of coalescence. Following a delay Tdeiay ~ 
T2 is exponentially distributed, as is apparent from the inset 
showing the cumulative distribution P(T2 > T). An exponen- 
tial e~~ 2TD / v is indicated as a black dashed line. Both plots 
contain simulation data for N = 10 4 (red), 10 5 (blue) and 
10 6 (green) as well as s = 0.01 (solid), s = 0.001 (dashed), 
and s — 0.0001 (dotted). For each parameter combination, 
random pairs are sampled at 10000 time points 2s~ 2 ^ 3 gener- 
ations apart. 



Assuming a molecular clock, the expected number of 
neutral differences between two genomes is ir = 2T2/i„, 
where fi n is the neutral mutation rate and T 2 is the time 
to the most recent ancestor of the pair of sequences. 
Across many realizations of the process (e.g. indepen- 
dent loci), T 2 follows a distribution P(T 2 ), which in the 
neutral case is exponential with mean N . Simulation re- 
sults of our model shown in Fig. [3] display a very different 
distribution of T 2 and equivalently it. Very few sequences 
coalesce early, which results in the long terminal branches 
observed in trees. We then observe a peak in coalescence 



around t 



2D 



lescence times decays exponentially with a characteristic 
time constant proportional to . Within a neutral coa- 
lescent frame work, a distribution of this kind would be 
interpreted as a rapid population expansion starting ^75 
in the past. Prior to this expansion, the population size 
would be estimated to have been constant at N e oc jjj. 
However, the size of the population did not change in our 
model. Instead, the population was adapting by many 
small steps, and the conclusion that N was decreasing in 
the past is wrong. 

Two lineages chosen at random from the population 
are most likely from near the center of the fitness dis- 
tribution. There are many individuals in this part of 
the distribution, so the probability of immediate coales- 
cence is therefore low. While the sampled individuals are 
typical, their ancestors tend to have higher than average 
fitness. Only after ancestral lineages have moved to the 
high fitness tail of the distribution, where only few indi- 



viduals are, does the rate of coalescence become appre- 
ciable. This migration of lineages towards higher fitness 
is a well kn own effect (|Hermisson et al.\ |2002[ |Q'fallon 



et al. 



2010| Rouzine and Coffin| |2007p and illustrated in 



after which the distribution of pair coa- i.e., T< 



Fig. [2 The speed at which lineages move toward higher 
(relative) fitness is initially the rate at which the mean 
fitness increases, while they slow down as they reach the 
tip. Consistent with the above interpretation, the delay 
of coalescence, Tdeiay , is roughly twice the time required 
for the mean fitness to catch up with the high fitness nose, 



delay 



2a-,. 



2D 



After lineages have moved to the 
high fitness tail, the seem to coalesce uniformly with a 
time constant T c ss 275. From the dependence of v on 
population parameters, we see that T c increases weakly 
with the population size as log 3 N where N — ND 3 . 



B. Site frequency spectra 

The density j(y)dv of neutral alleles in the frequency 
interval [v, v-\-dv\ is known as the site frequency spectrum 
(SFS). The SFS is a convenient summary of the neutral 
diversity segregating in the population. A mutation that 
happened on a particular branch of the genealogy will 
later be present in all individuals that descend from this 
branch. Hence the SFS harbors information about the 
distribution of branch weights and the branch length of 
the genealogy. In a neutral coalescent model, the site fre- 
quency spectrum is simply given by f(v) = where 
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© = is the average heterozygosity. Importantly, it 

is a monotonically decreasing function of the frequency. 
Fig. [4] shows site frequency spectra measured in simula- 
tions of our model. The most striking qualitative dif- 
ference is the non-monotonicity a feature known to be 



common in the presence of selective sweeps ( Fay and Wu 



2000 ) and in this context attributed to hitch-hiking. Sim- 



ilar phenomena arise in facultatively sexual organisms 



(Neher and Shraiman 2011). 



The non-monotonicity of f(v) implies the existence 
of long branches deep in the tree that are ancestors of 
almost everybody in the population, whereas a small 
minority of the population descends from different lin- 
eages. Such very asymmetric branchings are unlikely in 
a neutral coalescent model, where at any split the frac- 
tion of individuals that go left or right is uniformly dis- 
tributed (Wakeley 2008). Such asymmetric branchings 



are the rule in our model and frequently observed in re- 
constructed genealogies from rapidly adapting organisms, 
see Fig. [T] 

The axes in Fig. [4] are scaled to facilitate the com- 
parison to analytic results. At low frequenci es, the site 
frequency spectrum is proporti onal to v~ 2 (Basdevant 
and Goldschmidtl |2008[ |Neher and Shraiman 2011) and 



therefore much steeper than the neutral SFS f(v) oc v~ 
Hence even large samples will be dominated by single- 
tons. For high frequency derived alleles, we observe an 
increase of f(y) that is well described by 



/O) 



T c p, 



(i/-l)log(l-i/) 



(2) 



as v approaches one. 

The majority of the contributions to the increase of 
f(v) for v — ^ 1 stem from the very last coalescent event. 
In this last coalescent event, two or more lineages are 
merging. One of these lineages is typically the ances- 
tor of almost the entire sample, while the others share 
the remaining minority. The distribution of the offspring 
of these lineages and their number has been studied 



by G oldschmidt and Martin (Goldschmidt and Martin 
2005), who showed that the distribution of the size of 



the biggest lineage is asymptotically oc (1 — v)^ 1 . In the 
appendix, we derive the additional logarithmic prefactor. 
To compare the SFS of our model to that of the BSC 
across the entire range of v, we simulated the idealized 
BSC; comp. dashed line in Fig. [4j 

The non-monotonicity f{v) is a clear indication that 
the genealogies in this model with selection are funda- 
mentally different from neutral genealogies. In the King- 
man coalescent, neither constant or exponentially grow- 
ing population sizes give rise to non-monotonic SFS; see 
supplementary Fig. [8j 



C. The time to the MRCA 

In the neutral coalescent, the expected time to the 
most recent common ancestor of a sample of size n in- 
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FIG. 4 Site frequency spectrum of (derived) neutral alleles in 
rapidly adapting populations is non-monotonic with peaks at 
low frequencies and near fixation. The asymptotic behavior of 
the SFS at low and high derived frequencies is shown is shown 
as thick black lines. The dash-dotted black line is the allele 
frequency spectrum of the Bolthausen-Sznitman coalescent 
with N = 100000 averaged over 10000 runs (shifted downward 
by a factor of 10). 



creases only very slowly as the sample size is increased. 
This is a consequence of the even branching ratios of 
the neutral coalescent; an additional individual will most 
likely coalesce with existing samples before coalescing 
with the common ancestor of the existing sample. In 
particular, the average time to the most recent common 
ancestor of the entire population, Tmrca, is just twice 
the average pair coalescent time (assuming a constant 
population size) (Wakeley 2008). 

In simulations of our model of a rapidly adapting pop- 
ulation, we find that the average Tmrca increases slowly 
with N relative to the average pair coalescence time 
as shown in Fig. |7| T his is a generic property of the 
BSC (Berestycki |2009 ). The reason for this increase of 



Tmrca compared to T 2 is the very asymmetric branching 
of deep nodes in the BSC. Any pair of sequences is most 
likely a descendant of the same majority branch, which 
dominates the average T%. When sampling more indi- 
viduals, however, one is more and more likely to sample 
minority branches, which increases Tmrca of the sample 
relative to T 2 . 



III. ANALYSIS 

The simulation results presented above show that ge- 
nealogies arising in our model are clearly distinct from 
those expected in Kingman's coalescent and display 
a number of features reminiscent of the Bolthausen- 
Sznitman coalescent. We will now describe how this coa- 
lescent process emerges from the dynamics of the model. 

Individuals in our model have a heritable htness which 
determines the distribution of the number of immediate 
offspring they have. While fit individuals have on average 
more offspring than less fit ones, the fitness differences in 
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the population are small and the offspring distribution 
across the population is narrow. However, fitness is her- 
itable and fit individuals can have a very large number 
of distant great'-grandchildren. Hence the distribution of 
offspring after t generations, P(n, t) will be dominated by 
fit individuals and can have a very long tail. Eventually, 
there is exactly one individual that is the ancestor to the 
entire population, while all others leave zero offspring. 

As t increases, P(n, t) changes slowly from the initial 
narrow distribution to a broad distribution with a power- 
law at intermediate times; it then converges to N~ 1 5 n ,N 
at long times. The broad distribution at intermediate 
times is at the heart of the correspondence of genealogies 
in models of adapting populations and the Bolthausen- 
Sznitman coalescent. 

The BSC assumes that all individuals are exchangeable 
and that in every coalescent event a randomly chosen 
group of lineages merges into one. Each possible merger 
event has a specific rate associated with it, and the rate 
at which p individu als merge into one common ancestor is 
q p = T 2 _1 (p— (Berestycki 2009[ ) (the general expres- 
sion for the rates is given below in Eq. ([6])). In contrast, 
in the neutral coalescent, higher order coalescence is very 
rare oc TV - ^ -1 ). We will present the basic properties of 
the BSC briefly in the discussion. 

To appreciate how these coalescence rates can emerge 
from a model with selection, consider the number of in- 
dividuals rii that descend from an individual i that lived 
t generations in the past, with the additional constraint 
that n.j w N. The probability that p individuals sam- 
pled randomly from the population have a common an- 
cestor t in the past is then given by 



Q P (t) 




(3) 



where the average (.) is over all rij. Assuming that t is 
small enough that the different are still approximately 
independent, we can express Q p as 



Q P (t) 



N 
W)jo 



dzz p - 1 (e- zn ) N - 1 (n p e- zn ) 
dzz p ~ 1 e- N ' s ' At) (-l) p d p z ^ z (t) , 



(4) 



where we introduced the Laplace transform 1 — $ z (i) = 
S„ e~ zn P(n, t) and assumed N&%(t) < 1. In the ap- 
pendix, we show that Q(t) ~ z^m for t > Tdeiay = 
For a limited interval after t > Tdeiay, the probability 
that p individuals have a common ancestor increases with 
rate 



2D 1 
v p — 1 



(5) 



per unit time. Prior to Tdeiay, the rate of coalescence 
is very low. This is in agreement with Fig. [SJ where we 
found that little coalescence happened early, while co- 
alescence times are exponentially distributed after that 



with characteristic time T c v/2D for t > Tdeiay The 
relative rates of mergers of 2,3,. . . is consistent with the 
Bolthauscn-Sznitman coalescent, explaining our observa- 
tions for the frequency spectrum and the time to the most 
recent common ancestor. 

The branching process approximation used to derive 
the result Eq. ^ is valid only for short times but nev- 
ertheless gives us the relative rates of multiple mergers 
once coalescence sets in. For subsequent deeper coales- 
cent events, the relevant lineages are already at the tip 
of the fitness distribution, and this process repeats itself 
without the delay Tdeiay In fact, after this delay all re- 
maining lineages are in a narrow region at the tip of the 
fitness distribution shown in Fig. [5]4. The situation now 
resembles that of coalescence in FKPP waves: The fitness 
of the lineages are roughly equal, but they have to stay 
ahead of a moving front. We can therefore employ the 
phenomenological theory of genealogies in FKPP waves 

which confirms the above re- 
see Supplement. We 
mod- 



by (Brunet et al. 
suit 



2006) 



for the coalescent time scale; 
present an additional argument based on "tuned'' 
els introduced in ( |Hallatschek 2011 ) in the Supplement. 

To corroborate our analysis, we performed additional 
simulations that allow us to measure the Laplace trans- 
form of the distribution of pair coalescent times for very 
large populations. This Laplace transform with z = fi 
equals the probability of identity by state of two ran- 
domly chosen individuals. The algorithm is similar in 



2007), see sup- 



spirit to the algorithm by (Brunet et al. 
plementary information. 

Strictly speaking, the analogy to an exchangeable co- 
alescent model like the Bolthausen-Sznitman coalescent 
requires that different coalescence events one lineage un- 
dergoes are independent. For this to be true, individuals 
descending from a lineage have to distribute evenly across 
the fitness distribution c(x) between coalescence events, 
which requires a time T eq T c . This equilibration of fluc- 
tuations c an be explicitly studied in m odels of purifying 
selection (Neher and Shraiman 2012). The relaxation 
time of the distribution is roughly T c , and we should 
not expect a clean convergence against the Bolthausen- 
Sznitman coalescent. Nevertheless, we find it to be a very 
good model for the observed genealogies after account- 
ing for the delay. The underlying reason is that local 
equilibration in the region where the ancestral lineages 
are is fast (t « D- 1 / 3 ). This region, however, undergoes 
fluctuations on the time scale T eq , which modulate the 
overall rate of coalescence but don't affect the local dy- 
namics dramatically. The situation is more favorable for 
waves of FKPP type describing the spread of individuals 



in space, where T eq -C T c in large populations (Brunet 



et al. 2007). 



IV. DISCUSSION 

We have shown that in a simple model of adapting 
populations, the observed genealogies are qualitatively 
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FIG. 5 Panel A: The common ancestor of the population is 
confined to a narrow region in the tip of the fitness distribu- 
tion. The most recent common ancestor, i.e., a coalescence 
event, tends to be at even higher fitness. Panel B: The delay 
Tdeiay and the coalescence time scale T c can be extracted from 
the Laplace transform of P{T2) for very large populations and 
confirms that T c ~ ^j. 



different from the standard neutral coalescent. Instead, 
genealogical trees are characterized by long terminal 
branches and almost simultaneous coalescence of multiple 
lineages. At branching events deep in the tree, one com- 
monly observes that almost all individuals of the popula- 
tion descend from one branch, whereas very few descent 
from the other branches. Such skewed branching is un- 
likely in neutral coalescent models, regardless of the his- 
tory of the effective population size. One consequence of 
these uneven branching ratios is a non-monotonic site fre- 
quency spectrum (SFS) of derived neutral alleles. Com- 
pared to the neutral coalescent, the low frequency part 
of the SFS is much steeper, whereas the high frequency 
part shows a characteristic up-turn; see Fig. |4j These 
features are incompatible with the neutral coalescent. 

A given pair of lineages is unlikely to coalesce in the 
bulk of the fitness distribution. Typically, both lineages 
move into the high fitness tip of the population distribu- 
tion before they coalesce as illustrated in Fig. [2j This re- 
sults in long terminal branches and a distribution of het- 
erozygosities peaked at intermediate values. After this 
delay, the typical time to coalescence is again on the or- 
der of the time it takes the fittest individuals to dominate 
the population; see Fig. [3] In panmictic populations, this 
time depends on the logarithm of the population size and 
in our model is proportional to log 3 N . 

We argue that the exponential amplification of fit lin- 
eages is responsible for these observations and that coa- 
lescence in such rapidly adapting populations is generi- 
cally described by a modified Bolthausen-Sznitman coa- 



man 



lescent (BSC) (Berestycki 2009 



1998). The BSC is a specia 



Bolthausen and Sznit- 



caseof the large class 



of A-coalescent processes (Pitman 19991. Given the dis- 



tribution p(f) of the fraction / of the population that 
descends from a single individual in the previous genera- 
tion, the rate at which p out of k lineages merge is given 
by 



(6) 



A felP = / df P {f)F{\ - fy 



The Bolthausen-Sznitman coalescent corresponds to 
p(f) ~ f~ 2 for large /, in which case Eq. (|6| reduces 
to 



_ 1 (p-2)\(k-p)\ 
h - p T c (fc-1)! 



(7) 



or Eq. ([5| for the special case k = p. The total rate at 
which coalescence events happen is therefore 



k 



(8) 



in contrast to the neutral coalescent, where Xk oc k(k— 1). 
A coalescence event reduces the number of surviving lin- 
eages on average by log A: such that the average rate at 
which the number of lineages decreases is T~ 1 k\ogk. 
The average time needed to reach the common ancestor 
of a sample of size n is ~ T c loglogn, in contrast to 2T C 
in the Kingman coalescent. For a more in depth discus- 



sion, see the recent reviews by ( Berestycki 2009 Brunet 



and Derrida 2012). The BSC is easily implemented as 
a computer simulation by drawing an exponentially dis- 
tributed random number with mean A^ 1 to determine 
the time of the next event. The type of event is then 
chosen with probabilities proportional to Xk, P - 

The compatibility of a sample with the neutral coa- 
lescent model is typically asse ssed using statistics such 
as Tajima's D (Tajima 1989). Tajima's D compares 
the average number of pairwisc differences to the to- 
tal number of segregating sites in the sample. In the 
case of the Bolthausen-Sznitman coalescent, the average 
pairwise diversity is proportional to T2, while the total 
number of segregating sites is proportional to T271/ log 71 
(compared to T2 log n for the neutral coalescent) . This 
tremendous excess of segregating sites is a consequence 
of the very steep SFS at small frequencies and results 
in D oc — log(n). Similar results are expected for the F 

which 



1993) 



and D statistics by Fu and Li (Fu and Li 
compare mutations on internal and external branches. 

The models of adaptation we have studied have a nar- 
row offspring distribution. Nevertheless, the exponen- 
tial amplification of fit genotypes over many generations 
gives rise to a distribution of clone sizes with the required 
asymptotic behavior. The important lineages are those 
that run ahead of the distribution, expand faster, and 



take over a significant fraction of the population! Brunet 



et ah, 2006). Over even longer times, the fitness of an- 
cestors and descendants decorrelates, which allows one 
to approximate the genealogies with the abstract BSC, 
which assumes that there are no correlations in offspring 
number across generations. 

Conventionally, an increased variance in offspring num- 
ber is accounted for by defining an effective population 
size. With a clone size distribution p(f) ~ f~ 2 , however, 



the variance diverges with the population size ( Neher and 



Shraiman 



2011). Similar effects arise in other models 



with very skewed offspring distribution ( Eldon and Wake- 



ley 2006). As a consequence, the genealogies are domi- 



nated by rare anomalously large clones and described by 
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the BSC rather than the Kingman coalescent. The rate 
of coalescence is not set by TV -1 , but the rate at which 
clones expand and collapse. We would like to stress that 
evolutionary dynamics thereby remains highly stochas- 
tic, even in very large populations. Analogous behavior 
has recently been observed in models of individuals in- 
vading uninhabited territory (FKPP type waves) ( |Brunet | 
\et al.\ |2007|), ensembl es of super-critical branching pro- 



influx of deleterious mutations. Similarly, populations 
in a steady balance between deleterious and beneficial 



cesses ( Schweinsber g 2003) , recu rrent selective sweeps 
(Durrett and Schweinsberg |2005 ) , and models of facul - 



tatively sexual populations (Neher and Shraiman 2011 1. 



The BSC is not only a good model for genealogies of 
adapting populations, but also applies to populations un- 
der purifying selection in which Muller's ratchet clicks of- 
ten. The standard model for the distortion of genealogies 
by purifying selection assumes that deleterious mutations 
are rapidly purged and coalescence is neutral in the mu- 
tation free class with a reduced population size Ne~ u / S , 
where u is the deleterious mutation rate and s is the ef- 



fect size of deleterious mutations ( Charlesworth et al. 



1993). More elaborate analysis based on a fitness-class 



coalescent explicitly tracks lineages through the popu- 
lation and calculates the contribution to coalescence be- 



fore lineages reach the mutation free class( Walczak et al. 



mutations ( Goyal et al. 2012 ) have genealogies as found 
here for rapidly adapting populations. The reason for 
this qualitative difference is the fact that the nose of the 
wave is not steady, but constantly turning over. Differ- 
ent lineages are struggling to get ahead of everybody else, 
and, in the frame of reference of the population (that is, 
relative to mean fitness), are exponentially amplified. In 
contrast, dynamics of lineages in the mutation free class 
is neutral if Muller's ratchet does not operate. 

In the supplementary material, we show that the ar- 
gument that gave rise to the particular coalescence rates 
m Eq. (5} can be extended to a large class of models 
that are controlled by a small and fluctuating popula- 
tion of highly fit individuals. We thus argue that the 
Bolthausen-Sznitman coalescent generically emerges as 
a consequence of the exponential amplification of the 
clones corresponding to these highly fit individuals, to- 
gether with the seeding of novel lineages. The latter could 
happen via lucky diffusion to high fitness (our model 
here), via large effect beneficial mutations, or via lucky 
outcrossing. After some time, the distribution of lin- 
eage size follows a powerlaw with an exponent close to 



2012). However, all of this only holds as long as the muta- 2 (Desai and Fisher 2007 Neher and Shraiman 2011 



tion free class is large and Muller's ratchet does not oper- 
ate, which requires Nse~ u / S 3> 1 (Neher and Shraiman 



Yule 



this s 



1925). Given an effective offspring distribution of 



2012| |Stephan et qq|1993 |. 

Fig. [6j shows the SFS of derived neutral alleles for dif- 
ferent ratios u/s. For small u/s, the SFS are similar to 
those of the neutral coalescent with a reduced time to 
coalescence, in accordance with the background selection 
theory. However, as soon as the ratchet starts to click 
frequently, the SFS develop the non-monotonicity char- 
acteristic of the Bolthausen-Sznitman coalescent. 



rape, the Bolthausen-Sznitman coalescent follows 



(Brunet et al. 2007 Schweinsberg 2003). Whether the 



BSC also describes genealogies in scenarios where fitness 
is increased in rather large increments (compared to the 
population diversity) ( Gerrish and Lenski 1998 Schiffels 



et al. 2011 ) remains an interesting topic for future work 



Given the apparent universality of the Bolthausen- 
Sznitman coalescent in spatially expanding populations 
and panmictic adapting populations, it should be in- 
cluded as a prior in popular population genetic and phy 
logenetic inference programs such as BEAST (Drum- 




mond and Rambaut 2007) 



0.01 0.5 0.99 

derived allele frequency v 



FIG. 6 Site frequency spectrum of derived neutral alleles in 
a background selection scenario. As the ratio u/s is varied 
while keeping a = us constant, the SFS interpolate between 
the expectation for the Kingman and the BSC. 



If the ratchet is clicking fast, the population distri- 
bution resembles that of traveling wave models where 
selection on fitness variation can not keep up with the 
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Appendix A: Supplementary Information 



1. Frequency spectrum of derived alleles in the BSC 



A mutation occurring on some branch of the genealogy at time r will be present in all individuals that descend from 
the branch considered, i.e., it's copy number equals the size of the T-family associated with that branch. These family 
sizes follow a Poisson-Dirichlet distribution PD(0, e~ T ) and the size z of the first block is distributed as (Berestycki 
[2009] ) 



P(z,t) 



r(i-e- T )r( e - r ) 



(Al) 



9 



Common alleles, i.e., those present in the majority of the population, most likely fall into this block, 
averaged frequency spectrum of those alleles is therefore approximately given by 



The time 



/(*) 



r(l -e- T )T{e- T ) 



(A2) 



1 z-y(\-z)y- x 

y yT(l-y)T(y) 
1 sin(7r#)z _ »(l - z^ 1 

where we have used the Euler reflection formula for the Gamma functions. The term sin(7ry)/7rj/ is close to unity for 
the small y where the dominant contribution to the integral comes from. Hence we have approximately 



The integral 



/(*)« / d yz -y(i-z)y- 1 

JQ 



2z- 1 



z(l-z)(logz + log(l-z)) 



(A3) 



dy 



sin(7ry)z v (l 



\y-i 



(A4) 



has an expression in terms of special functions, and higher order corrections to Eq. (A3) can be obtained by expanding 
the sine. They do not change the result qualitatively. 

This expression is expected to be the dominant contribution to the frequency spectrum of common alleles, at 
1 — z -C 1. In this limit, it describes the simulation data very well, see Fig. [4j At the rare end of t he spectrum , 
the observed power law f(z) ~ 



is a (well-known) consequence of exponentially growing populations (Yule 19251, 



corresponding to expanding clones in the high-fit tip of the traveling wave. 



2. Generating function of the offspring number distribution 

The generating function of the distribution of the number of descendants of individuals initially at fitness x is 
defined as <j) z {x, t) = 1 — ^2 n P(n, t, x)e~ nz . The initial condition P(n, 0, x) — S nt % translates into <p z (x, 0) = 1 — e~ z . 
The generating function obeys the following equation 



d t cf> z (x,t) = Dd 2 x 4> z (x,t) - vd x <j> z {x,t) + x<f> z (x,t) - <^(x,t) 



(A5) 



It is easy to see that for small z this equation has two simple asymptotic solutions with a sharp crossover separating 
the two. 



(f> z (x,t) 



X < X 

x > x 



(A6) 



The cross-over is approximately at x* = y — ^ lQg Z J X . Note that the branching process description is only 

valid for a limit time before population size constraints have to be imposed and lineages cease to be approximately 
independent (Cohen et al. 2005 Hallatschek 2011). At short times, however, it is a valid approximate description 

vt 2 

of the number of offspring of a marked individual. The mean, for example, is given by d z <p z (x,t) e xt a" at short 
times. We are interested foremost in the distribution of offspring averaged over the population c(x). For short times 
with x c < x* we have 



3>(z,t) = z Ae "2 



dxAi 



. ( v 2 



x 



(A7) 



We see that the behavior of this average changes qualitatively as t > ^ when it starts to be boundary dominated. 
Assuming x* < x c , i.e., the saturation of (f> z (x,t) happens within the validity of the Airy function solution, we have 



$(z,i) fa zAe~ 



dxAi 



( v 2 



\4D3 



x 



c (*-2Tj) 



+ x*A 



dxAi 



AD 3 



x 

Di 



(A8) 
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Both of these integrals are dominated by a narrow region in the vicinity of x* . Plugging in the definition of x* , we 
find 



$(z,()«i*4e 2D 



dSxAi 



( v 2 



Sx 



Di 



e '*(*-A) + / dSxA 



Sx 



4Z?s 



e 2r> 



= C(t,a;*)e"^ = C{t,x*)e 2Df +^-— = E{t,x*)z^ 
where C(i, x*) and E{t,x*) depend only weakly on x*. The pth derivative of this is 

t) « *-*$(*,*) JJ(-^- -t) = (-!)**-*$(*, t) J! (i 



i=0 



_v_, 
2Dt' 



We can now evaluate Qp. 



Q P (t) 



N 

m Jo 



«te« p_1 e- v *^(-l) 1 '^(« ) t) 



2£>t ' 



r(p) 

2^nCo X (*'-5B*) 



dzz- 1 e~ Ar$(2) 7V$(z,i) 



nr=i( 



Choosing i = fg, we have 



r(p) 



>(*) 



2Dt i 



r(p) 



r(p - r" 



r(i - r-i)r(p) 



(A9) 



(A10) 



(All) 



(A12) 



For large r, all rates become equal, which simply means that the population has coalesced into its common ancestors. 



For t = + St, we can differentiate in St to 



q P = d t Q p {t)\t= v / 



2D 



2D 1 

v 1 — p 



(A13) 



which are the Bolthausen-Sznitman merger rates. 

We have assumed a steady distribution with a well-defined cutoff at x c . The actual high fitness end of the population 
fluctuates quite a bit, which would result in large population size fluctuations, if the population size was not tightly 
controlled via a mean fitness that is catching up. This feed back from the mean fitness, however, is delayed and does 
not affect the dynamics at the nose instantaneously. The mapping to the FKPP models discussed below is better 
suited to describe the effects of fluctuations of the nose on coalescence. 



3. Relation to FKPP models 



After the equilibration time T eq the ancestors of all extant individuals are located in a region of width oc D 1 / 3 in 

2 

the nose of the fitness distribution roughly x c rs ahead of the mean. Within this region, the growth rate does not 
vary much and the deterministic population distribution is a solution of 



D 



dc(x) 
dx 2 



9c(x) 

v— 1- x c c{x) = 

ox 



(A14) 



which has an exponential solution c(x) ~ e X1 with v 
corrected by a cut-off ( Brunet and Derrida 19971) , which corresponds to v 



( Cohen et al. 2005 1 . Brunet et al 



is of Bolthausen-Sznitman type wit 



~ + Dj. Waves of this from select the minimal velocity 



2\/x c D and 70 = 2D/v, consistent with 
Brunet et al.\ |2006| |2007[) have argued that the coalescence process of such waves 



1 a typical time scale given by 



T, = 



log n 

7T 2 7^"( 7o ) 



21og J 



(A15) 
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where n is the size of the population relevant for the FKPP wave. In our case, only the part of the wave that has a 
substantial chance of being the common ancestor of the population is relevant and this part is located in an area of 
with 0(D^ 3 ) away from the tip. The number of individuals in this range is 

r° _^ v 

log noc logy 1/3 dxe 2 ° 21)2/3 ( Al6 ) 

Substituting this into the above and using x c = v 2 /4D, we find 

T 2 oc £ (A17) 

in agreement with the above findings. Furthermore, this argument explains why we continue to see Bolthausen- 
Sznitman coalescence after the branching process description is no longer valid. 



4. The BSC in tuned models 



The calculation giving rise to Eq. ^ was done for a specific model of selection and mutation. We now want 
to argue that similar considerations hold for a large class of related models. The evolution and composition of 
traveling evolutionary waves is governed by the fitness distribution in the population, c(x), and the probability that 
an individual at fitness x will take over the entire population u(x) (Hallatschek 2011[ ). The deterministic and linear 
parts (i.e., up the a cutoff) of u{x) and c(x) are left and right eigenfunctions of the same evolution operator C The 
distribution c(x) is a localized pulse, while u(x) is an increasing function of x (high fitness individuals have a higher 
chance of spreading than low fitness individuals). 

Looking forward in time, u(x) is the probability that a lineage at x is the common ancestor of the entire population 
in the distant future. Before the descendants of a single lineage take over the population, they spread evenly across 
the fitness distribution. Once spread out, the chance of fixation of this lineage is given by the fraction of off-spring it 
contributes to the total population. Hence a fraction u(x) of the population will typically trace back to an individual 
at x a time T eq in the past, where T eq is the equilibration time of individuals across the population (T eq v/2D in 
our case). 

After identifying u(x) with the fraction of individuals descending from a single ancestor at x, the probability that 
p randomly chosen individuals have a common ancestor T eq generations in the past is simply 



Qp (^"et 



dxc(x)u p (x) 



(A18) 



Since u(x) is monotonic in x we can change variables to u 



(T eq ) fa / du c{u) 



du\ 
dx 



(A19) 



If the evolution operator C is dominated by selection, u(x) grows faster than exponential in the relevant range of x 
and furthermore is proportional to the inverse of c(x). Hence Eq. (A20) reduces to 



Qp (T e , 



du u 2 u p 



1 



(A20) 



The constant of proportio nality depends on t he details of the stochastic dynamics which we have calculated above 
by other means (see also (Hallatschek 2011) for an alternative calculation of Q2(T eq ) — J dxcu 2 ). According to 



this reasoning, the fitness of common ancestors of the population are distributed as c{x)u{x). Using this relation to 
compute u{x) from the observed distribution of ancestors, we can predict the distribution of most recent common 
ancestors that were the root of a 2-merger event: their fitness should be distributed as c(x)u 2 (x), which compares 
well with simulations. 



5. Laplace transform of pair coalescence times 



Instead of simulating individuals, one can also simulate a population spread out on a 2 dimensional grid, where 
one dimension corresponds to fitness and the other is a neutral direction used to study the pairwise diversity of 
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individuals. This keeps a record of the pair-coalescence time distribution: The distribution of pairwise differences 
between individuals can be decomposed into the distribution of coalescence times and distribution of distances along 
the neutral axis, given the coalescence time T. The latter follows a diffusion equation with diffusion constant D = /i/2. 
Hence the averaged distribution over all pairwise differences along the neutral coordinate should be 



Hence the Fourier transform of the pair-wise distance function is the Laplace transform of the pair-coalescence time 
in the variable k 2 D/2. 

We know from the individual based simulations that the coalescence process is essentially broken into two phases, 
a waiting period where nothing is happening, followed by an exponential decay. The Laplace transform of such a 
function is 




(A21) 



Fourier transforming this distribution results in 




(A22) 



P(z) = 



l + T c z 



(A23) 



These two parameters can be easily extracted by fitting this function to simulation results. 
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FIG. 7 The time to the most recent common ancestor increases relative to the average pair coalescent time with the population 
size N = ND 1 ' 3 , a well known feature of the Bolthausen-Sznitman coalescent. 
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FIG. 8 The allele frequency of neutral derived mutations in the Kingman coalescent with constant (a = 0) and exponentially 
growing population size (a > is the growth rate measured in units of N^ 1 ). The genealogies are produced with the program 
ms ([Hudson 2002]). The black lines show the theoretical expectation for a population of constant size and a rapidly expanding 
population. The x-axis is scaled as in Fig. [4] 



